function   res = momentsClosedForm(model,setupFilter)

% For efficient starting values when computing unconditional second moments
global Var_z2old Var_z3old
if ~exist('Var_z2old','var')
    Var_z2old = [];
end
if ~exist('Var_z3old','var')
    Var_z3old = [];
end
% Construct model structure to match the model output
modelGMM = doModelforGMM(model);

% The moments in the normal distribution
ne         = size(model.eta,2);       %number of shocks
vectorMom3 = zeros(1,ne);
vectorMom4 = 3*ones(1,ne);
vectorMom5 = zeros(1,ne);
vectorMom6 = 15*ones(1,ne);

% Computing moments
sig         = 1;
maxAuto_lag = setupFilter.CSregress_m;
mom2nd      = UnconditionalMoments_2nd_Lyap(modelGMM.gx,modelGMM.gxx,modelGMM.gss,...
    modelGMM.hx,modelGMM.hxx,modelGMM.hss,modelGMM.eta,sig,vectorMom3,vectorMom4,maxAuto_lag,Var_z2old);
Var_z2old = mom2nd.Var_z;
if setupFilter.appOrder == 3
    mom3rd = UnconditionalMoments_3rd_Lyap(modelGMM.gx,modelGMM.gxx,modelGMM.gss,...
        modelGMM.gxxx,modelGMM.gssx,modelGMM.gsss,...
        modelGMM.hx,modelGMM.hxx,modelGMM.hss,modelGMM.hxxx,modelGMM.hssx,modelGMM.hsss,modelGMM.eta,sig,...
        vectorMom3,vectorMom4,vectorMom5,vectorMom6,mom2nd.Var_xs,mom2nd.Var_xfxf,mom2nd.Var_xs_xfxf,maxAuto_lag,Var_z3old);
    Var_z3old = mom3rd.Var_z;
end
if setupFilter.appOrder == 1 || setupFilter.appOrder == 2
    % hours and wages are in deviation from steady when matched to the data
    pos_l      = find(strcmp(modelGMM.labely,'$l_t$'));
    pos_w      = find(strcmp(modelGMM.labely,'$w_t$'));
    mom2nd.E_y(pos_l,1) = mom2nd.E_y(pos_l,1)-modelGMM.g0(pos_l,1);
    mom2nd.E_y(pos_w,1) = mom2nd.E_y(pos_w,1)-modelGMM.g0(pos_w,1);
    
    % recall Var(y) = E[y^2]-E[y]E[y'], so E[y^2] = Var(y)+E[y]E[y']
    tmp               = diag(mom2nd.Var_y+(mom2nd.E_y+modelGMM.g0)*(mom2nd.E_y+modelGMM.g0)');
    EyEy2             = [mom2nd.E_y(1:setupFilter.nyMom,1) + modelGMM.g0(1:setupFilter.nyMom,1);...
                        tmp(1:setupFilter.nyMom)];
    
    % Loadings for yields in the extended state space system (with pruning)
    D_yields         = 4*[model.byx  model.byx  model.Byxx];
    
    % Loadings for the slope in the CS regressions in the extended state space system (with pruning)
    D_slopeCS        =  D_yields(setupFilter.CSregress_m+1:end,:)-D_yields(setupFilter.CSregress_m,:);
    
    % Getting Cov(yields^k,slope^k)
    Cov_y_k_slope_k = D_yields(setupFilter.CSregress_m+1:end,:)*mom2nd.Var_z*D_slopeCS';

    % Getting Cov(yields^(k-m),slope^k). Recall Cov_z(:,:,m) = Cov(z_t+m,z_t)
    Cov_y_kminusM_slope_k  = D_yields(1:end-setupFilter.CSregress_m,:)*...
                             mom2nd.Cov_z(:,:,setupFilter.CSregress_m)*D_slopeCS';

    % Getting Var(slope^k) 
    Var_slope              = D_slopeCS*mom2nd.Var_z*D_slopeCS';                         
    
    % First and second moments of all yields
    meanYields             = 4*model.by0 + D_yields*mom2nd.E_z+4*(0.5*model.byss*sig^2);
    stdYields              = sqrt(D_yields*mom2nd.Var_z*D_yields');
    
    % First and second moments of shrinkage moments used in estimation
    meanForEstim           = modelGMM.g0+mom2nd.E_y;
    stdForEstim            = sqrt(diag(mom2nd.Var_y));
elseif setupFilter.appOrder == 3
    % hours and wages are in deviation from steady when matched to the data
    pos_l      = find(strcmp(modelGMM.labely,'$l_t$'));
    pos_w      = find(strcmp(modelGMM.labely,'$w_t$'));
    mom3rd.E_y(pos_l,1) = mom3rd.E_y(pos_l,1)-modelGMM.g0(pos_l,1);
    mom3rd.E_y(pos_w,1) = mom3rd.E_y(pos_w,1)-modelGMM.g0(pos_w,1);

    % recall Var(y) = E[y^2]-E[y]E[y'], so E[y^2] = Var(y)+E[y]E[y']
    tmp               = diag(mom3rd.Var_y+(mom3rd.E_y+modelGMM.g0)*(mom3rd.E_y+modelGMM.g0)');
    EyEy2             = [mom3rd.E_y(1:setupFilter.nyMom,1) + modelGMM.g0(1:setupFilter.nyMom,1);...
                        tmp(1:setupFilter.nyMom,1)];
    
    % Loadings for yields in the extended state space system (with pruning)
    D_yields         = 4*[model.byx+3/6*model.byssx*sig^2 model.byx model.Byxx model.byx 2*model.Byxx model.Byxxx];
    
    % Loadings for the slope in the CS regressions in the extended state space system (with pruning)
    D_slopeCS        =  D_yields(setupFilter.CSregress_m+1:end,:)-D_yields(setupFilter.CSregress_m,:);
    
    % Getting Cov(yields^k,slope^k)
    Cov_y_k_slope_k = D_yields(setupFilter.CSregress_m+1:end,:)*mom3rd.Var_z*D_slopeCS';
    
    % Getting Cov(yields^(k-m),slope^k). Recall Cov_z(:,:,m) = Cov(z_t+m,z_t)
    Cov_y_kminusM_slope_k  = D_yields(1:end-setupFilter.CSregress_m,:)*...
                             mom3rd.Cov_z(:,:,setupFilter.CSregress_m)*D_slopeCS';

    % Getting Var(slope^k) 
    Var_slope              = D_slopeCS*mom3rd.Var_z*D_slopeCS';
    
    % First and second moments of all yields
    meanYields             = 4*model.by0 + D_yields*mom3rd.E_z+4*(0.5*model.byss*sig^2 + 1/6*model.bysss*sig^3);
    
    stdYields              = sqrt(diag(D_yields*mom3rd.Var_z*D_yields'));
    
    % First and second moments of shrinkage moments used in estimation
    % hours and wages are in deviation from steady when matched to the data
    meanForEstim          = modelGMM.g0+mom3rd.E_y;
    stdForEstim           = sqrt(diag(mom3rd.Var_y));
end

Cov_excessReturn_slope = Cov_y_kminusM_slope_k-Cov_y_k_slope_k;
idx                    = setupFilter.CSregress_m+1:size(D_yields,1);
scaling                = setupFilter.CSregress_m./(idx-setupFilter.CSregress_m);
CS_cov                 = [NaN(setupFilter.CSregress_m,1);scaling'.*diag(Cov_excessReturn_slope)];
CS_loadings            = [NaN(setupFilter.CSregress_m,1);diag(Cov_excessReturn_slope)./(scaling'.*diag(Var_slope))];
%CS_cov                 = [NaN(setupFilter.CSregress_m,1);diag(Cov_excessReturn_slope)];
%CS_loadings            = [NaN(setupFilter.CSregress_m,1);diag(Cov_excessReturn_slope)./(scaling'.*diag(Var_slope))];
%% Output
if setupFilter.SMMwithCSreg == 2
    % Including CS moments and adjusted CS moments
    %risk-adj CS moments most be the empirical moments and therefore added later
    CS_var           = [NaN(setupFilter.CSregress_m,1);scaling'.^2.*diag(Var_slope)];
    res.modelMoments = [EyEy2; CS_cov(setupFilter.maturites_CS);...
                        CS_var(setupFilter.maturites_CS)]';
elseif setupFilter.SMMwithCSreg == 1
    %res.modelMoments = [EyEy2; CS_cov(setupFilter.maturites_CS)]';
    CS_var           = [NaN(setupFilter.CSregress_m,1);scaling'.^2.*diag(Var_slope)];
    res.modelMoments = [EyEy2; CS_cov(setupFilter.maturites_CS);CS_var(setupFilter.maturites_CS)]';
else
    res.modelMoments = EyEy2';
end
% Removing first moment for mean of hours and wages
res.modelMoments    = res.modelMoments(3:end);

res.CS_cov           = CS_cov;
res.CS_loadings      = CS_loadings(setupFilter.CSregress_m:setupFilter.CSregress_m:end)';
res.meanYields       = meanYields;
res.stdYields        = stdYields;
res.meanForEstim     = meanForEstim;
res.stdForEstim      = stdForEstim;
scalingAll           = [NaN(1,setupFilter.CSregress_m) scaling];
res.scaling          = scalingAll(setupFilter.maturites_CS);


end